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Abstract 

In this paper we explore new ways to study the zero temperature limit of quantum statistical 
mechanics using Quantum Monte Carlo simulations. We develop a Quantum Monte Carlo method 
in which one fixes the ground state energy as a parameter. The Hamiltonians we consider are of 
the form H = Hq + XV with ground state energy E. For fixed Hq and V, one can view as a 
function of A whereas we view A as a function of E. We fix E and define a path integral Quantum 
Monte Carlo method in which a path makes no reference to the times (discrete or continuous) at 
which transitions occur between states. For fixed E we can determine \{E) and other ground state 
properties of H. 
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I. INTRODUCTION 



Quantum Monte Carlo methods are widely used to compute properties of quantum sys- 
tems using classical sampling algorithms. In this paper we develop a novel Quantum Monte 
Carlo method that allows one to numerically investigate ground state properties of a quan- 
tum system. 

A virtue of Quantum Monte Carlo is that one is not required to manipulate vectors in the 
Hilbert space corresponding to the quantum system. The dimension of this Hilbert space 
typically grows exponentially with the physical size of the system. Instead, Quantum Monte 
Carlo methods map the problem of approximating the ground state energy (or some other 
observable) onto the problem of evaluating an expectation value with respect to a probability 
distribution q{X) over a set of configurations C (so X G C). In order to evaluate this 
expectation value, one can use a classical Markov chain Monte Carlo algorithm to sample 
configurations from the distribution q. Markov chain Monte Carlo works by defining a 
Markov chain on the space of configurations C. This Markov chain can be described by 
an update rule which tells you how to generate a new configuration of the chain from the 
current one. The Markov chain is constructed so that the limiting distribution is q{X). 
One then applies some large number A^o of iterations of the Markov chain to some initial 
configuration Xq. If Nq is sufficiently large then after these iterations, the distribution of 
subsequent configurations will be arbitrarily close to q. 

We now give a brief description of how our method is used to estimate properties of the 
ground state. We write the Hamiltonian as H{X) = Hq + XV^ where Hq is diagonal in a 
given basis {\z)} and XV is off diagonal in this basis. (In an n spin system 2; is an n bit 
string.) In section [ll] we outline our assumptions and restrictions on Hq and V. With these 
choices, the ground state energy is always less than or equal to zero, and we will see that 
for each value of < 0, there exists one positive value X{E) such that the ground state of 
H{X{E)) has energy E. To use our Monte Carlo method, one first must fix E < and a 
large integer m. We define a path of length m to be a sequence {zi, ...,Zm}, where each n 
bit string Zi is the label of the state \zi). These paths are the configurations of the previous 
paragraph. We will define a probability distribution f{{zi, Zm}) over the set of all paths 
of length m (this distribution is also a function of the value of E which was chosen). We 
will show how the function X{E) can be obtained by computing an average with respect to 
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the probability distribution /. 

To motivate our method and to get a general idea of how it works, consider the function 



G{E,\) =Tr 



-X 



V 



(1) 



H{X)-E 

Assuming that E < and A > are chosen so that the Taylor series expansion converges, 
we can write 

-A \ 1 



G{E,X) = Tr 



1 + 



Hn-E 



V Ho-E 



V 



m=l 



-A 



:V 



(2) 



Hq — E 

It is clear from the expression in equation [T] that the function G{E, A) blows up when 
E Eg{X), where Eg{X) is the ground state energy of H{X). Equivalently we can say 
that at a fixed value of E the blow up occurs as A ^ ^{E)., where Eg{X{E)) = E. At this 
value of A the Taylor series expansion must diverge. In fact, this divergence occurs because 
as m becomes large, terms in the series approach 1 for large m (here we have made some 
assumptions about the Hamiltonian which we discuss in the next section) so 

-HE) 



Tr 



Hn-E 



V 
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For the remainder of this section we assume that m is large enough to make ~ close to =. 
By inserting complete sets of states in the basis that diagonalizes Hq we can express the 
LHS as a sum over paths 

m ^ 

{X{E)r Yl {^l\-V\^m){z„.\-V\Zm^,)...{z2\-V\z^)l[^—^^l 

{zi,...,Zrn.} « = 1 * 

where Ei = {zi\HQ\zi). Now taking the log and differentiating with respect to E, we obtain 
1 dX{E) 



{zi,...,Zm} 



^'\m^E,-E 
^1=1 ' 



X{E) dE 

i=l 

Here the expectation value is taken with respect to the measure / on paths defined by 

1 m ^ 

f {{Zi, ...,Zm}) = —{Zi\ - V\z^){Zm\ - V\Zm-l)...{z2\ - V\zi) JJ 



1=1 



where F is a normalizing constant. We will show how to sample with respect to the distri- 
bution / in a way that makes numerical work possible. Sampling from the distribution / 
will also allow us to compute — x^^^^ from equation jsj as well as other properties of the 
ground state. 

Our paper is organized as follows. In section [IT] we describe the types of Hamiltonians 



for which our method apphes. In section III we outhne the new method that we propose. 



In section [IV] we explicitly construct Monte Carlo update rules for the case where V = 
— J2'i=i ^'^d we give numerical data using our algorithm at n = 16 where we are able 
to compare with exact diagonalization. In section |V] we review the continuous imaginary 
time Quantum Monte Carlo method [6], which is based on the thermal path integral. We 
also derive a novel estimator in this ensemble of paths for the ground state energy which 
becomes exact in the limit (3 —>■ oo. 



II. THE HAMILTONIAN 



We consider finite dimensional Hamiltonians of the form 

H{\) =Ho + XV, 

where Hq is diagonal in a given basis {\z)}, and V has zeros along the diagonal in this basis. 
We make the following assumptions about the Hamiltonian: 

1. The off diagonal matrix elements of V in the basis {\z)} which diagonalizes Hq are 
all either negative or zero. (This ensures that our Quantum Monte Carlo method will 
not suffer from a sign problem.) 

2. The ground state of H{X) is not degenerate for any value of A G {—oo, oo). 

3. The smallest eigenvalue of Hq is zero. Note that this condition can be fulfilled without 
loss of generality by adding a constant term to the Hamiltonian. Writing \zo) G {\z)} 
for the unique state with Ho\zo) = 0, we further require that V\zo) ^ . (This implies 
that \zq) is not an eigenvector of V since (-ZolV'l-Zo) = follows from assumption 1 
above.) 
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We write I'lpgiX)) and Eg{X) for the ground state eigenvector and ground state energy of 
H{X) . From second order perturbation theory in A, we have that 



^9 



dX^ 





V\zo)\' 




He 





< (4) 



dX 



where the inequahty is strict because V\zo) 7^ 0. 
Using the fact that 

dE 

= {MX)\V\4jg{X)) (5) 

we show that 

f 

> , for A<0 
= , for A = 
< , for A>0 . 

In order to obtain the inequahties, we use the variational principle. When A > 0, the 
ground state energy must be less than zero, since \zq) has zero expectation value for H 
(and is not an eigenvector of H{X)). This, together with the fact that Hq is positive 
semidefinite, imphes that {ipglVlijjg) < 0. The analogous result for A < is obtained in the 
same way. These inequalities give a qualitative picture of the curve Eg{X). Starting from 
Eg{0) = 0, the curve slopes downwards as it goes out from A=0, and approaches —00 on 
both sides of the origin for sufficiently large |A|. Note that this imphes that for each E <0 
there is one positive and one negative value of A (call them X{E) and X^{E) respectively) 



such that Eg{X{E)) = E and Eg{X-{E)) = E. Furthermore, we show in appendix Al that 
it is always the case that 

X{E) < \X^{E)\ . (6) 

We refer to the case where the inequality is strict as the generic case. We illustrate the 
qualitative features of the curve Eg{X) (for the generic case) in figure [T] In the nongeneric 
case where equality holds at some particular value of E , then in fact equality holds at every 
value of E and the curve -^^(A) is symmetric about A = 0. 
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X 




Figure 1: Eg{X) for the Hamiltonians we consider. As A ^ ±00 we have Eg —>■ —00. 
III. A NEW QUANTUM MONTE CARLO METHOD 

Definition of the Ensemble of Paths and Relevant Estimators 

As motivated in the Introduction, we now define an ensemble where the configurations 
are sequences {zi, ...,Zm} (where each Zi is an n bit string) , and we show how properties of 
the ground state can be computed in this ensemble. We refer to the sequences {zi, ...jZm} 
as paths. 

To begin, we fix < and a large integer m as parameters. As in section|ll| we take X{E) 
to be the positive value of A such that H{X) has ground state energy with corresponding 
eigenvector \iljg{\{E))). We now describe how our method allows us to approximate X{E) 
and other properties of the ground state. 

Recall from the Introduction that the probability distribution / over paths is defined by 

1 m ^ 

'""^ ^ F{E,m) ^''^ - ^1^-) • • • - V\z,) n ^-r^ (7) 
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with^ 



F{E,m)= Yl {^i\-V\^m)...{z2\-V\z^)l[ — 



{zi,...,Zm} 



Tr 



i=l 



E,-E 



Hn — E 



vy 



(8) 
(9) 



As examples, we now define two quantities ^'^^J^'^ and X'^{E,m) as ensemble averages with 
respect to the distribution / on paths 



/3{E,m) = ^ f{{zu...,Zm})f3est{{zi,...,Zm}) = {Pest)f 

{zi,...,Zm} 

X^{E,m) = f{{zi,...,Zm})Xlsti{Zl,--,Zm}) = {\lst)f 

{zi,...,Zm} 

where we have defined the estimators (hence the subscript) 



(10) 



m ^ 

Pesti{zi, Zm}) = 



i=l 



E,. - E 



1=1 ^ ' 



(11) 

(12) 



with Ei = {zi\Ho\zi), and Zm+i = zi, Zm+2 = Z2- Our reason for using the symbol Pest will 
become clear in section |V] where we will discuss its interpretation as an inverse temperature. 
We show in appendix A that the ensemble averages ^'^^^^ and X^{E,m) correspond to 
properties of the quantum ground state in the hmit m ^ oo 

1 d\{E) 



m-»oo m 



\{E) dE 
1 1 



_ \{E){^lj,{E)\V\^,{E)) 
lim A2(E,m) = {\{E)f . 



(13) 

(14) 
(15) 



Equation 13 is equation |3] of the introduction and equation 14 follows from [5j One can also 
derive expressions for higher derivatives of log(A(i?)) as averages with respect to /. 



^ The nongeneric case where equaUty holds in equation |6]can arise when there exists a unitary transformation 
U such that U'^VU = —V and WHqU — Hq. In this case it is seen from equation [o] that F{E,m) = 
when m is odd. In the nongeneric case m must always be taken to be even. 
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Monte Carlo Simulation 



We propose to use the measure / as the basis for Monte Carlo simulations. In particular, 
our Monte Carlo algorithm begins by choosing E < and a large integer m and then 
samples sequences of bit strings from the distribution /. One can then use the estimators 



from equations 10 and 12 to evaluate the quantities X{E) and {iljg{X{E))\V\iljg{X{E))), using 



equations [14] and [15] for the limiting behaviour of these estimators. We do not construct a 
general method for sampling from /, instead we leave it to the reader to construct such a 
method for the particular choice of V at hand. We do note that it is essential that whatever 
method is used conserves the total number m of transitions in the path. We now give an 
example of this for a generic spin system with V = — (y\. 



IV. AN EXAMPLE: TRANSVERSE FIELD SPIN HAMILTONIANS 

We now explicitly construct a method to compute averages with respect to the distribu- 
tion / in the case where the Hilbert space is that of n spin \ particles, and V = — 'YTi=\ 
The Hamiltonian i^o is an arbitrary diagonal matrix in the Pauli z basis for n spins. Our 
algorithm can be used to compute the average of any function of the path which is invari- 
ant under cychc permutations of the path. For this choice of V ^ the spectrum of -ff (A) is 
symmetric about A = (so this corresponds to the nongeneric case where equality holds in 
equation ]6] for all E < 0). With these choices, the paths which have nonzero weight (with 
respect to /) are periodic paths of m bit strings of length n where each string z differs from 
the previous one by a bit flip. We must take m to be even since each bit must flip an even 
number of times so that the path is periodic (and so the total number of bit flips m in the 
path must be even). In order to sample from these paths according to /, we construct a 
Markov Chain which has / as its limiting distribution (actually, our Markov Chain con- 
verges to the correct distribution over equivalence classes of paths which are only deflned 
up to cyclic permutation-this is why we restrict ourselves to estimating quantities which are 
cychcally invariant). Note that with our choice of V we can write (see equation [T]) 

1 " 1 

Our Markov chain is defined by the following update rule which describes how the configu- 
ration is changed at each step 
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1. Choose an integer i G {1, ...,m} uniformly at random. 

2. Consider the bit strings Zi-i, Zi, Zi-^i in the current path (where Zm+i = zi). Suppose 
that Zi-i and Zi differ in bit qi G {1, ...,^}, which we write as Zi = Zi^i © e^^. Also 
write q2 G {1, ...,n} for the bit in which Zi and Zj+i differ, so Zj+i = -Zj © e^a- 



3. If gi 7^ g2 then propose to change the bit string Zi to the new value Zi = Zi-i 
Accept this proposal with probability 



E, - E 



Precept = min 1 1, ^ 

where Ei = (zj|ifo|^j)- This Monte Carlo move has the effect of interchanging 2 
consecutive flips in the path (see flgure[2]). 

4. If qi = q2 , then choose a new bit qnew ^ {1, ■■.,n} from the probability distribution 

Piqnew = q) = ^^^ (16) 

where E' = {zi^i © eq\Ho\zi-i © Cg), and W = X^^i e'^-e - Then (with probability 
1) change Zi to the new value Zi_i © e^^^^. This Monte Carlo move replaces a pair of 
consecutive flips which occur in the same bit with 2 new flips in a possibly different 
bit (see figure |3]). 

We show in appendix [C] that this algorithm can be used to estimate any quantity which is 
invariant under cyclic permutations of the path (note that all estimators we have discussed 
have this property). 



Numerical Simulation with a Particular Choice of Hq 

We have numerically tested our new Monte Carlo algorithm using a C++ computer 
program. In this section we show numerical data at 16 bits where we are able to compare 
results with exact numerical diagonahzation. We studied the Hamiltonian with V as in the 
previous section, and Hq corresponding to the combinatorial optimization problem Exact 
Cover which stems from our interest in quantum computation [3] 

1=1 
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Figure 2: Monte Carlo update where the order of 2 flips in the path is interchanged 



i+1 



i+1 



-i-l" 



Figure 3: Monte Carlo update where 2 adjacent flips in the path which occur in the same bit are 
replaced by flips in a different bit. 



where 



E 

c=l 



cr 



ii(c) 



42 (c) 



Here Hq is a sum over Nc terms, which in computer science are called clauses. Each clause 
involves three distinct bits zi(c), Z2(c), ^3(0). A clause is said to be satisfied by an n bit 
string z if the state \z) has zero energy for the corresponding term in the Hamiltonian. The 
particular choice of and the bits involved in each clause defines an instance of Exact 
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Cover. Such an instance is said to be satisfiable if there is an n bit string Zg which satisfies 
all clauses in the instance. In that case the z basis state l^^) is the zero energy ground state 
of i7o, that is Hq\zs) = 0. 

We generated an instance of Exact Cover on 16 bits with a unique satisfying assignment 
through a random procedure. Figures [i] and [5] show the values of \{E) and —\{E)dE{X) / d\ 
computed using equation [TO] for 200 values of E, with m = 1000. Statistical errors were 
computed using Ulh Wollf's error analysis program [10]. This data set was taken by running 
10* Monte Carlo updates on each of two processors of a dual core laptop computer for each 
value of E. The two processors ran simultaneously and the total time taken for all the 
data was under 5 hours. We also include the curves for these quantities obtained by exact 
diagonalization, which are in good agreement with the Monte Carlo data. Since it is hard 
to see the error bars in figures |4] and [5| we have plotted the errors separately in figures [6] 
and[7l 
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Monte Carlo Data and Exact Diagonalization for X{E) 
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Figure 4: X{E) computed using Monte Carlo data and exact diagonalization for a 16 spin Hamil- 
tonian. Statistical error bars are included for the Monte Carlo results, but they are barely visible. 
We have also plotted the errors separately in figure |6] 
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Monte Carlo Data and Exact Diagonalization for <-'k V > 

18 I \ 1 1 \ 1 \ 1 I "T" 



Exact Diagonalization 




Figure 5: —X{E){ipg{X{E))\V\Tpg{X{E))) computed using Monte Carlo data and exact diagonaliza- 
tion for a 16 spin Hamiltonian. Statistical error bars are included for the Monte Carlo results, but 
they are barely visible. Errors are also plotted separately in figure [7| 
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statistical Error Estimates for X(E) 



o ^*V%i(\o 
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o Actual Error 

+ Estimated Statistical Error 




Figure 6: A(-E') from figure [i} The black crosses show the estimated statistical error. The blue 
circles show the magnitude of the difference between the Monte Carlo estimates and the result of 
exact numerical diagonalization. 
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Figure 7: —X{E){^lJg{X{E))\V\'^pg{X{E))) from figure |5j The black crosses show the estimated 
statistical error. The blue circles show the magnitude of the difference between the Monte Carlo 
estimates and the result of exact numerical diagonalization. 
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V. A NEW ESTIMATOR FOR THE ENERGY IN STANDARD CONTINUOUS 
IMAGINARY TIME QUANTUM MONTE CARLO 



In the ensemble of paths which was considered in the previous section, the parameters 
m and E are fixed and then quantities l3{E,m) and X^{E,m) are computed as ensemble 
averages. In this section we review the standard approach to computing thermal averages 
using continuous imaginary time Quantum Monte Carlo [6j. To use this method, parameters 
/? and A are fixed beforehand and quantities E{(3, A) and m(/3, A) are computed as ensemble 



averages. We will also derive a new estimator for {H) 



Tr[He 



Z((3) 



which is valid for large 



/3 simulations. The form of this estimator establishes a connection between the continuous 
imaginary time Quantum Monte Carlo method and the new method that we outlined in the 
previous section. 

The standard method |6] is based on the expansion of the partition function 



Tr 



Tr 



Tr 



-I3H 



^(_A)-e-^^o / dt^ / dtm-l... / ^ 
^-n ^0 Jo Jo 



dtiVl{trr,)Vl{tm^l)...Vl(ti 



m=0 



m=l 



/3 



dtjYi I dtfYi_i. 



Xr Yl {^l\VM{Zm\V\z^.l)...{z2\V\Zi) 

{zi,...,Zm} 

t2 1 
^^^^-{Eltl+E2(t2-tl) + ...+El{l3~tr^)) 



(17) 

'0 ^0 Jo 

where Viif) = e^Hoy^-tHo ^nd Ei = {zi\Ho\zi). This expression is interpreted path 
integral, where a path is defined by a piecewise constant function z{t) for t G [0,/?]. The 
function z{t) takes values in the set {z} which are the labels of the basis states {1^)} which 
diagonahze Hq. In the above expression for Z{P), we have 



z{t) 



Z2, 



0<t <ti 

ti<t<t2 
t-m—l ^ t ^ tr, 



The allowed values of m are 2, 3, 4, 5, ... as well as m = in which case the path is constant 
z{t) = Zi (for some Zi) for all t G [0,/3]. 
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We define Ho{z{t)) to be {z{t)\Ho\z{t)) , so Ho{z{t))dt = Eih + ^2(^2 - ^i) + ••• + 



Ei{(3 — tm) ■ Equation 17 can be used to define a measure p on paths in imaginary time. 



The probability of a given path parameterized by z{t) is 



1 „-|3Hn{z^) ?77 — fl 

Z{P)^ , III — V 

^A™(zi| - V\z^){zJ - V\z^.i)...{z2\ - V\zi)e~So^o(z{mt^^^_^^^^ ^ _^ 

where Z{P) is the normalization. Our assumption that off diagonal elements of V are 
nonpositive guarantees that p > 0. The task of samphng paths from the distribution p can be 
accomplished using Markov chain Monte Carlo methods such as those outlined in references 
OmiH]. By using these methods to sample paths from this probability distribution, one can 
compute physical properties of the quantum system at an inverse temperature /3. Known 
estimators for the expectation of the terms in the Hamiltonian are (see appendix [B]) 



(^0) = ^ ={^J H,{z{t))dt), (18) 

Here m is the number of transitions in the path z{t) and is not fixed. Note that the estimator 
for {W) only involves the number of transitions in the path. The above two estimators can 
be combined to obtain (H). One can also obtain the following expressions for the variances 
of these estimators in the limit /3 — > 00 
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This shows that as /3 — > cxd, the distributions for these estimators become sharply peaked 
about their mean values. We define intensive ensemble averages 

rh{P,X) _ m 

E{P,X) ^ [' Ho{z{t))dt),-^^^ 



which in the — 00 limit become respectively the ground state expectation value of —XV 
and the ground state energy. 
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A New Estimator for the Ground State Energy 



We now derive the following novel estimator for the energy in the standard ensemble, 
which is useful in the hmit /3 — > cxd : 



{H) = {E^), + O 



where E* is a function of the path defined to be the smallest value of E which satisfies the 
equation 

m+l ^ 



i=l 



In this expression the Ei are the energies of the states \zi) visited along the path, with 



Em+i = El. Note that the above equation is almost identical to equation 11 (this justifies 
our choice of notation Pest)- We obtain this formula by a similar method to that used in 
reference |T] to obtain an alternate estimator for the energy (H). Our formula, however, is 
valid when some or all of the {E^} are the same, and therefore resolves a serious difficulty 
encountered in reference p^. This may be of use in large P Monte Carlo simulations, as an 
alternative to the standard estimator for {H). We note in particular that this estimator 
does not involve the times {ti, ...,tm} in the path. 

Equation [20] can be derived by considering the Laplace transform of the partition function 
(with s > -Eg) 



H + s 



/ e-^'Z{P)dp = Tr 
Jo 

We can also express this as 

/•oo /"OO °° r 

/ e-^'Z{P)dp = / dPe-^'T E {^i\-yM{zm\-V\z^_i)...{z2\-V\zi) 

Jo Jo L r_ _ 1 



m=0 L {zi,...,Zm} 
fOO 



/»oo fOO m+l 

Jo Jo 



m+l 



Ei + s 



= E^'" E i^^\-V\^m)...{z2\-V\z,)l[^ 

m=0 {zi,...,Zrri} i=l 

Performing the inverse Laplace transform gives 

oo -, n / m+l 



oo If/ "'"T^ 1 

m = j:>^"' E i-^\-y\-m)...{z.\-v\z,)— dse^^in — 

m=0 {z^,...,z^} ''^ ^ i=l ' 
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C is a contour in the complex plane which encircles the poles of the integrand, which are 
located at {—Ei}. The expectation value of the energy can then be expressed as 

1 



Tr[He 



-0H] 



1 dZ 
Zdi3 



/3s 



ds 



2Tvi JC ^ \ lli=l 



rm+1 1 



ds 



(21) 



The complex function h{s) = e 



m+l 1 



Ei+s 



will in general have multiple saddle points 



along the real axis. We can solve for the locations of these saddle points by writing 



h{s) 



Ms) 



where 



m+l 



s — 



I J2 log(^. + s) 



i=l 



Saddle points occur at real values s* where ^(s*) = , which says that 



m+l 



1 



In the integrals in equation [2T| we can choose the contour of integration along the curve 
of steepest descent through that saddle point s* which is largest (it is possible to show that 
one can deform the contour to follow this curve without changing the encircled poles). 

Performing the integrals in equation [21] and letting E* = —s* we obtain 



where E* is the smallest solution to equation 20 



VI. CONCLUSIONS 



We have outlined a new approach to Quantum Monte Carlo simulations in which prop- 
erties of the ground state of a quantum system are computed at a fixed value of the ground 
state energy. We have confirmed the validity of our method by performing a numerical 
simulation of a system consisting of 16 spins. Our approach involves a path integral which 
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does not include any jump time variables, as in the stochastic series expansion [71 18]. We 
have also obtained a new estimator for the ground state energy which is valid in continuous 
imaginary time Quantum Monte Carlo simulations but which does not involve the imaginary 
time variables. We hope that our method will be applied to numerical simulations that go 
beyond the toy system studied in this paper. 
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Appendix A: DERIVATION OF ESTIMATORS FOR THE NEW MONTE CARLO 
METHOD 

1. Properties of the operator A{E) 

Our analysis of the new Monte Carlo method is based on properties of the operator 



We take E < Q, and \{E) is defined to be the positive value of A such that the ground 
state of H{\) has energy E. We will show the following properties of this operator: 

1. All eigenvalues oi A{E) are real and < 1 in absolute value. 

2. \'4>g{\{E))) is an eigenvector of A{E) with eigenvalue +1. It may also be the case 
that \'^g{—\{E))) is an eigenvector of A{E) with eigenvalue —1. There are no other 
eigenvectors of A{E) with eigenvalues ±1. 

The second property says that the subspace of states spanned by eigenvectors of A{E) with 
±1 eigenvalues is either 1 or 2 dimensional. We will see that the case where it is two 
dimensional occurs only when H{—X{E)) has ground state energy E. 

To show that eigenvalues of A{E) are real whenever E < 0, note that the spectrum of 
A{E) is the same as the spectrum of the operator 



1992. 




as was seen from the Introduction. 



^/lh^A{E) 



1 



^JHq — E 
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which is Hermitian. 

Now suppose (to reach a contradiction) that |r) is an eigenvector of A{E) with eigenvalue 
R>1. Then 

r 1 1 

|r) = E\r) 



Ho + \{E)^V 



which says that there exists an eigenvector of if (^^) with eigenvalue E. In section II we 
proved that the ground state energy of H{\) is strictly decreasing for positive A, which 
means that no eigenvector of if (^^) can have energy smaller than or equal to E. This is 
a contradiction and so all positive eigenvalues of A{E) are < 1. 

Now, since we have proven that all eigenvalues of A{E) are < 1, if it is the case that 
some negative eigenvalue of A{E) is < — 1 then the eigenvalue W of A{E) which is largest 
in magnitude has negative sign. If this were true, then the limit 



lim Tr 



2k+l 



would be equal to a negative constant. But this cannot be the case since all matrix elements 
of A{E) are positive or zero. So we have shown property 1. 

We now proceed to show property 2 which was stated earlier in the appendix, and in the 
process we prove the inequality |6] from section [llj The real nonzero eigenvalues oi A{E) can 
be related to eigenvalues of ii(A) for some value of A. In particular, suppose that a; is a real 
nonzero eigenvalue of A{E) with eigenvector \uj). Then 

{Ho-E)uj\uj) = -X{E)V\uj) 

so 

Let us write the eigenvalues of A{E) as 

1 = ai{E) > a2{E) > ... > a2n{E) > -1 . 

(which follows by property 1.) Then the values of A at which ii(A) has an eigenvalue with 
energy E are 

m 

aj{E) 
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for j G {1, ...,2"}. Write X-{E) for the negative value of lambda such that H{\^{E)) has 
ground state energy E. The values of A which are smallest in magnitude (and hence closest 
to the axis A = 0) are -^j^ = \{E) > and -^^^ = ^-(E) < and these correspond to 
the ground state at energy with |A_(i?)| > A(-E'). All other values are greater than X{E) 
in magnitude. This proves inequality [6| and also shows property 2 described above. 



2. Derivation of the Estimators Pest and Ag^^ 

Having derived properties 1. and 2. of the operator A{E) in the previous section, we 
now proceed to use these properties to prove equations [14] and [15] 

Our treatment below apphes to both the generic case where inequality [6] is strict as well 
as the nongeneric case where equality holds as long as in the latter case m is always even. 
In either case we have 

Tr[{A{E)y'^] = {\{E)rE{E,m). (Al) 
(Recall the definition of F{E,m) from equation [s]) We can write this as 

2" 

ixiE)rF{E,m) = j2iamr ■ 
1=1 

Taking the log and differentiating both sides gives 



2"^ 

[log (A(E)-) + log {FiE, m))] = \ Yl ^ Mmr-' ^ • (A2) 

Now take the limit as m ^ cxd. Note that in the nongeneric case the fact that m is taken 
to be even ensures that the denominator of equation X2 does not vanish in the limit of 
large (even) m. The limit of the RHS is zero for any fixed value of E. This is because in 
the large m limit the only terms which contribute to the sum in the numerator are those 
corresponding to values of j for which \aj{E)\ = 1. For these values oiaj{E) (note j is either 
1 or 2" for these values) it is always the case that ^ = so these terms contribute zero to 
the sum. 

So we have shown that 

f I d\ 1 dF(E,m)\ , 

lim m— — — + — — \ ' M = 0. A3 

m-oo \^ \{E)dE F{E,m) dE J ^ ^ 
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Hence 



1 1 dF{E,m) 



lim 

m^oo m F{E, m) dE 



1 d\ 



\{E) dE' 

The left hand side of this equation can be rewritten as an ensemble average, which results 
in _ 

1 d\ 



lim 



dF{E,m) 



m. — >nci m 



m^oo 777, F(E, m) dE m^oo m ^ '-^'Ei-E'^ \{E)dE' 

So for fixed m sufficiently large, one can estimate the quantity — y7^44 using the ensemble 



(A4) 



\{E) dE 



14 



average ^^^^^ ■ This proves equation 

We now derive an estimator for the quantity \{E) itself as an ensemble average. For this 
purpose we make use of the fact that 

F{E,m- 2) 



lim 

m.-*oo F{E,mj 



x{Ey. 



Expanding the numerator and denominator as sums over paths, we obtain 

F{E,m-2) _ Ew,...,._,}(^i|-^k™-2)...(^2|-V^ki)nr=f ^ 



rm 1 



FiE,m) Ew,......}(^il - V\zJ...{z,\ - V\z,) UT=i e.-e 

Now rewrite the numerator as 



(A5) 



F{E,m-2) = ^ (^{Zl\-V\Zm){Zm\-V\Zm-l){Zm^l\-V\Zm-2)---{z2\-V\Zi) 



{zi,...,Zm.} 

m ^ 



i=l 



E,. - E 



Sz^zm-iiEm - E){Ei - E)- 



(^il^'ki). 

where 5z^zm-i is the Kronecker delta. Since our distribution f{zi,...,Zm) is invariant under 
cychc permutations of the bit strings {zi}, we can write 

F{E,m-2) = ^ i{zi\-V\z^){zm\-V\zm-i){zm~i\-V\zm-2)---{z2\-V\zi) 

{zi,...,Zm} 

m r m ll\ 



where Zm+i = z\ and 2:^+2 = ^2- Inserting this formula into equation X5 gives the final 
expression for A(i?)^ as an ensemble average 



lim 



m-+oo \ m 



ra 



{zi\V'^\zi) I J 



X{Ef 



(A6) 



This proves equation 15 
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Appendix B: ESTIMATORS IN CONTINUOUS IMAGINARY TIME QUANTUM 
MONTE CARLO 



In this section we derive the known estimators for (Hq) and {XV) stated in equations 18 
andHa 



Estimator for (Hq) 



To derive the estimator for {Hq) we write (using the Dyson series to expand e 



Tr[Hoe-^^] _ 1 



Tr 



m=0 



t2 



where Vj{t) = e^^^Ve (The m = term in the above sum is ^^^'"[e ^^°].) Inserting 
complete sets of states in the basis {1^)} which diagonahzes Hq we obtain 

Tr\H p~l^^] 1 °° r 

^i:r^W = zTM E (-^r E {zi\Ho\zi){z,\V\z^){zjV\zm-i)...{z2\V\z,) 



Jo JO 



dtr. 



dt 



t2 



^^^^-{Eiti+E2{t2-t{) + ...+Ei{l3-tm)) 



'Mz{t = mp- 



where the expectation value is with respect to the measure p defined in section |V| and 
Ho{z{t = 0)) = {z{0)\Ho\z{0)). Noting that the measure p is invariant under a translation of 
the path by a time x G [0, /3] (this corresponds to the transformation U goes to (U+x) mod /3 
for i G 1, m followed by a reordering of the labels i to maintain time ordering), we have 
that 

{Ho{z{t = 0))), = {Ho{z{t = x))), , for all xg[0,/5]. 
We obtain the stated estimator for {Hq) by averaging over all x G [0,(3] 

Tr[Hoe-^^] 



{Ho) 



1 

/, Ho{z{x))dx) p, . 
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Estimator for {XV) 



As in the previous section, we begin by expanding the operator e 



Tr[XVe-^^] _ 1 



Tr 



m=0 



Jo Jo Jo 

(Bl) 



t2 



For m = 0, 1, 2... we have 



I dtm+l / dtm--- / dtiVl{tm+l)Vl{tm)...Vl{tl)S{ti) 

'o Jo Jo 

dtm+l / dtm... / dt2Vi{tm+l)Vi{tm)--Vlit2)V 

10 Jo Jo 

f>/3 r-tm nt2 

dtm / dtm-i... / c?tiV7(t„)V}(t™_i)...V7(ti)V^. 

'0 Jo Jo 



Plugging this expression into |Bl| we obtain 

Tr[XVe-^"] 1 



Z{(3) 



Tr 



-l)Y^{-Xr^'e 



■/3Ho 



m=0 

tm + l 



t2 



dtm+l / dtm--- / dtiVj(tm+l )Vjitm)-Vjit,)5it^] 

Jo Jo 

oo 



Zi(3) 



Tr 



m=l 

rtm rt2 
dtm I dtm~l--- I dtiVi{tm)Vj{tm-l)---Vj{ti)6{ti) 

Jo Jo 

= -((l-5m,0)5(tl))p 

m 

= —((1 — 6m,o) ^ ^{ti))p (since only ti can ever be 0). 

1=1 

In the last two hues of the above m appears inside an expectation value {---)p- In this context 
m is considered to be a function of the path. Now we can use the fact that the measure p 
over paths is invariant under translations of the path in imaginary time to write 

TrlXVe-^"] 1 ^ 

j ^ = -((l-5,n,o)^ / Y.^{t,-t)dt), 

P Jo 

= -(-) 



Tr[e-P"] 



So {XV) is —4 times the average number of jumps in a path of length (3. 
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Appendix C: CONVERGENCE OF THE MARKOV CHAIN 



We show in this section that the Markov Chain defined in section [IV] can be used to 
estimate any quantity which is invariant under cyclic permutations of the path. In order to 
streamline the proof, it will be useful to define a different Markov Chain over paths which 
has the update rule (for some fixed < p < 1) 



1. With probability p do 1 update of the Markov Chain defined in section IV 



2. With probability 1 — p apply a random cychc permutation to the path by letting 
{zi,Z2,...,Zm} {zj,Zj+i, ...,Zm,zi, ...Zj_i} fov Uniformly random j G {l,...,m}. 

In the next two sections we show that the above Markov Chain has limiting distribution / 
(defined in equation [?]) for any choice of the parameter < p < 1. This Markov Chain with 
fixed < p < 1 induces a random walk on equivalence classes of paths where an equivalence 
class is the set of all paths related to a given path by cyclic permutation. In this equivalence 
class random walk, step 2 does nothing. So the hmiting distribution over equivalence classes 
is the same whether or not step 2 is performed. If one only estimates quantities which are 
invariant under cyclic permutations (note that the estimators we have discussed have this 
property) then one can use the algorithm with p = 1 (so step 2 is never performed). 

To show that the above Markov Chain converges to the limiting distribution / over paths, 
it is sufficient to verify that the update rule constructed above satisfies the following two 
conditions [5]: 

• Ergodicity: Given any two paths A and B, it is possible to reach path B by starting 
in path A and applying the Markov chain update rule a finite number of times. 

• Detailed Balance: For any two paths A and B, 

f{A)P{A -.B) = f{B)P{B ^ A) (CI) 

where P{X Y) is the probability of transitioning to the path Y given that you start 
in path X and apply one step of the Markov chain. 

We now show that this Markov Chain satisfies these conditions. 
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Ergodicity 



In order to show ergodicity of the Markov Chain defined above, we first note that a path 
can be specified either by a fist of bit strings {^i, ...,Zm}, or by one bit string Zgtart followed 
by a fist of bits in which fiips occur {61, bm}, with each br G {1, n}. 

We now show that by applying the Monte Carlo update rules illustrated in figures [2] 
and [sj it is possible to transform an arbitrary path A < — > {zstart, {bi, bm}} into another 
arbitrary path B < — > {VstartS^i, Cm}} where i/start and Zstart differ by an even number of 
bit fiips. This is sufficient to show ergodicity because any path B can be cychcally permuted 
into a path which starts in a state ystart that differs from Zgtart by an even number of fiips 
(and our Markov chain includes moves which cyclically permute the path). We assume here 
that m > 4, since we are interested in the limit of large m anyways. In order to transform 
path A into path B, we give the following prescription: 

1. First transform Zstart into y start- To do this, note that one can move any two fiips bi and 
bj so that they are just before and just after Zstart (i.e 61 = bj and bm = bi), by applying 
the fiip interchange rule illustrated in figure [2] So long as you do not interchange the 
first and last fiip in the fist, the bit string Zgtart will remain unchanged. Then one can 
fiip both of these bits in the bit string Zgtart by interchanging the two fiips 61 and bm- 
This describes how to fiip any two bits in Zstart, assuming that fiips in these bits occur 
somewhere in the path. Now suppose that you wish to fiip 2 bits in Zstart but one or 
both of the bits does not occur in the current list of fiips in the path. In that case you 
must first take some pair of fiips which occur in some other bit q, and then move them 
until they are adjacent using the fiip interchange rule (without ever moving them past 
Zstart)- Once they are adjacent, you can replace them with a pair of fiips in another 
bit using the fiip replacement rule illustrated in figure [3]. Assuming m > 4, there will 
always be two pairs of fiips in the path which can be replaced by fiips in the two bits 
that you desire to change in Zstart- 

2. After Zstart has been transformed into ystart, one must then make the list of fiips equal 
to {ci,...,Cm}. This can be done by interchanging fiips and replacing pairs of fiips as 
described above. Since this can always be achieved without interchanging the first and 
last fiip in the path, the bit string ystart will remain unchanged by this procedure. 
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Detailed Balance 



Here we demonstrate that the Markov Chain defined above satisfies the detailed balance 
condition from equation|Cl[ To show this, fix two paths A and B and consider the probability 



of transitioning between them in one step of the Monte Carlo update rule. This probability 
is zero except in the following cases 

1. A = B. In this case detailed balance is trivially satisfied. 

2. A cychc permutation of the bit strings (which is not the identity) maps the path A 
into the path B. In this case f{A) = f{B) and the probability of transitioning from A 
to B in one move of the Markov Chain is also equal to the probability of the reverse 
transition from B io A (this is because for every cyclic permutation which maps A 
to B the inverse permutation is also cychc and maps B io A). So detailed balance is 
satisfied. 

3. A and B are the same path except for at one location. In other words, A can 
be described by the sequence {zf^...,z:^} and B can be described by the sequence 

, z^} where the corresponding bit strings are all the same except for one pair 
zf and zf. Write qi for the bit in which zf and zfi-^ differ, and q2 for the bit in 
which zf and z^-^ differ. We have to consider two cases depending on whether or not 
qi = q2- 

• Case 1: gi 7^ g2- Detailed balance follows in this case since the transition probabilities 

follow the Metropolis Monte Carlo rule 

P(A^B) ( Ef-E^ 1 

^ ^ min U * ' 



fjA) ^ Ef-E _ PjB^A) 
f{B) Ef -E P{A^ B) 



Case 2: qi = q2. In this case, from equation [16] we have 

P{A-^B) _ Ef-E 
P{B A) ~ Ef-E 

and 

f{A) _ Ef-E _ PjB^A) 
f{B) Ef -E P{A^ B) 
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